home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / NLEqn.cc < prev    next >
C/C++ Source or Header  |  1996-03-03  |  4KB  |  197 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #if defined (__GNUG__)
  24. #pragma implementation
  25. #endif
  26.  
  27. #ifdef HAVE_CONFIG_H
  28. #include <config.h>
  29. #endif
  30.  
  31. #include "NLEqn.h"
  32. #include "dMatrix.h"
  33. #include "f77-fcn.h"
  34. #include "lo-error.h"
  35.  
  36. extern "C"
  37. {
  38.   int F77_FCN (hybrd1, HYBRD1) (int (*)(int*, double*, double*, int*),
  39.                 const int&, double*, double*,
  40.                 const double&, int&, double*,
  41.                 const int&);
  42.  
  43.   int F77_FCN (hybrj1, HYBRJ1) (int (*)(int*, double*, double*,
  44.                     double*, int*, int*),
  45.                 const int&, double*, double*, double*,
  46.                 const int&, const double&, int&,
  47.                 double*, const int&);
  48. }
  49.  
  50. static NLFunc::nonlinear_fcn user_fun;
  51. static NLFunc::jacobian_fcn user_jac;
  52.  
  53. // error handling
  54.  
  55. void
  56. NLEqn::error (const char* msg)
  57. {
  58.   (*current_liboctave_error_handler) ("fatal NLEqn error: %s", msg);
  59. }
  60.  
  61. // Other operations
  62.  
  63. int
  64. hybrd1_fcn (int *n, double *x, double *fvec, int *iflag)
  65. {
  66.   int nn = *n;
  67.   ColumnVector tmp_f (nn);
  68.   ColumnVector tmp_x (nn);
  69.  
  70.   for (int i = 0; i < nn; i++)
  71.     tmp_x.elem (i) = x[i];
  72.  
  73.   tmp_f = (*user_fun) (tmp_x);
  74.  
  75.   if (tmp_f.length () == 0)
  76.     *iflag = -1;
  77.   else
  78.     {
  79.       for (int i = 0; i < nn; i++)
  80.     fvec[i] = tmp_f.elem (i);
  81.     }
  82.  
  83.   return 0;
  84. }
  85.  
  86. int
  87. hybrj1_fcn (int *n, double *x, double *fvec, double *fjac,
  88.         int *ldfjac, int *iflag)
  89. {
  90.   int nn = *n;
  91.   ColumnVector tmp_x (nn);
  92.  
  93.   for (int i = 0; i < nn; i++)
  94.     tmp_x.elem (i) = x[i];
  95.  
  96.   int flag = *iflag;
  97.   if (flag == 1)
  98.     {
  99.       ColumnVector tmp_f (nn);
  100.  
  101.       tmp_f = (*user_fun) (tmp_x);
  102.  
  103.       if (tmp_f.length () == 0)
  104.     *iflag = -1;
  105.       else
  106.     {
  107.       for (int i = 0; i < nn; i++)
  108.         fvec[i] = tmp_f.elem (i);
  109.     }
  110.     }
  111.   else
  112.     {
  113.       Matrix tmp_fj (nn, nn);
  114.  
  115.       tmp_fj = (*user_jac) (tmp_x);
  116.  
  117.       if (tmp_fj.rows () == 0 || tmp_fj.columns () == 0)
  118.     *iflag = -1;
  119.       else
  120.     {
  121.       int ld = *ldfjac;
  122.       for (int j = 0; j < nn; j++)
  123.         for (int i = 0; i < nn; i++)
  124.           fjac[j*ld+i] = tmp_fj.elem (i, j);
  125.     }
  126.     }
  127.  
  128.   return 0;
  129. }
  130.  
  131. ColumnVector
  132. NLEqn::solve (int& info)
  133. {
  134.   ColumnVector retval;
  135.  
  136.   int n = x.capacity ();
  137.  
  138.   if (n == 0)
  139.     {
  140.       error ("equation set not initialized");
  141.       return retval;
  142.     }
  143.  
  144.   double tol = tolerance ();
  145.  
  146.   retval = x;
  147.   double *px = retval.fortran_vec ();
  148.  
  149.   user_fun = fun;
  150.   user_jac = jac;
  151.  
  152.   if (jac)
  153.     {
  154.       Array<double> fvec (n);
  155.       double *pfvec = fvec.fortran_vec ();
  156.  
  157.       int lwa = (n*(n+13))/2;
  158.       Array<double> wa (lwa);
  159.       double *pwa = wa.fortran_vec ();
  160.  
  161.       Array<double> fjac (n*n);
  162.       double *pfjac = fjac.fortran_vec ();
  163.  
  164.       F77_XFCN (hybrj1, HYBRJ1, (hybrj1_fcn, n, px, pfvec, pfjac, n,
  165.                  tol, info, pwa, lwa));
  166.  
  167.       if (f77_exception_encountered)
  168.     (*current_liboctave_error_handler) ("unrecoverable error in hybrj1");
  169.     }
  170.   else
  171.     {
  172.       Array<double> fvec (n);
  173.       double *pfvec = fvec.fortran_vec ();
  174.  
  175.       int lwa = (n*(3*n+13))/2;
  176.       Array<double> wa (lwa);
  177.       double *pwa = wa.fortran_vec ();
  178.  
  179.       F77_XFCN (hybrd1, HYBRD1, (hybrd1_fcn, n, px, pfvec, tol, info,
  180.                  pwa, lwa));
  181.  
  182.       if (f77_exception_encountered)
  183.     (*current_liboctave_error_handler) ("unrecoverable error in hybrd1");
  184.     }
  185.  
  186.   if (info < 0)
  187.     retval.resize (0);
  188.  
  189.   return retval;
  190. }
  191.  
  192. /*
  193. ;;; Local Variables: ***
  194. ;;; mode: C++ ***
  195. ;;; End: ***
  196. */
  197.